A Measurement of the Appearance of Atmospheric Tau Neutrinos by Super- Kamiokande 
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Super-Kamiokande atmospheric neutrino data were fit with an un-binned maximum likelihood method to 
search for the appearance of tau leptons resulting from the interactions of oscillation-generated tau neu- 
trinos in the detector. Relative to the expectation of unity, the tau normalization is found to be 1.42 ± 
0.35 (stat) Iq'J 4 ( s y s ) excluding the no-tau-appearance hypothesis, for which the normalization would be zero, 
at the 3.8<T level. We estimate that 180.1 ±44.3 (stat) ^js' 8 (sys) tau leptons were produced in the fiducial vol- 
ume of the detector by tau neutrinos during the running period, compared to an expectation of 120.2 (sys) 
events. 
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It is now well known that neutrinos undergo flavor oscil- 
lations. The flavor states of the neutrino measured through 
the weak interaction are quantum mechanical mixtures of neu- 
trino mass states. As observed in the quark sector, this mix- 
ture results in oscillations of detected flavor states. Evidence 
exists for this effect in atmospheric neutrinos 0] solar 
neutrinos EH7), reactor experiments @, and long baseline 
oscillation experiments lt9l4TTl. In 2011, the T2K |12|, MI- 
NOS 1 13 1, and Double Chooz [14| experiments showed the 
first indications of full three-flavor oscillations. In 2012 the 
Daya Bay lfT31 and RENO IfTBI experiments reported the first 
precision measurements of the 613 mixing angle which drives 
three-flavor oscillation. 

Definitive proof of flavor oscillation requires unambiguous 
appearance of the charged current interaction of a neutrino 
not in the original source. In the dominant oscillation for 
at GeV energies, — > V T oscillations, observing the result- 
ing T lepton is quite difficult. This is because producing a tau 
lepton requires a neutrino of energy greater than a threshold 
of 3.5 GeV. Long -baseline experiments tuned to the neutrino 
oscillation maximum for their distances tend to have the bulk 
of their neutrinos below this energy. Furthermore, the tau lep- 
ton immediately decays to final states with an electron, muon 
or mesons plus a tau neutrino so the tau lepton itself cannot be 
easily seen. Nevertheless, the OPERA collaboration was re- 
cently able to show evidence for a single reconstructed event 
in their emulsion consistent with tau appearance lfl7l , The 
Super-K collaboration first published a search for tau appear- 
ance in atmospheric neutrinos in 2006 fl8l . Since the atmo- 
spheric neutrino flux extends to energies well above 10 GeV, 
and spans a wide range of baselines, we expect to see tau 
leptons produced in the Super-K detector. However, these 
events must be distinguished from other high-energy atmo- 
spheric neutrino interactions. Further comparisons of these 
techniques can be found in lfT9"l and prospects for future de- 
tectors can be found in l20l . 

This letter reports a result from a new search utilizing 
the Super-Kamiokande experiment. Super-Kami okande (Su- 
per-K) is a 50,000 ton water Cherenkov detector[21| with 
22.5 ktons of fiducial volume. It consists of two concentric de- 
tectors: a inner-detector with 11,129 inward-looking 20 inch 
photo-detectors, and an outer-detector with 1885 outward- 
facing 8 inch photo-detectors which acts as a veto. Its large 
target mass makes it well suited to look for the rare appear- 
ance of tau neutrinos from oscillations. The typical energy of 
atmospheric neutrinos is about 1 GeV. Due to the previously 
noted energy threshold, our expectation is that about one v T 
charged current event per kton-year will be produced in the 
Super-K detector. 

Super-K has been in operation for approximately fifteen 
years and has had several running configurations indicated 
by the labeling SK-I (1996-2001), SK-II (2002-2005), SK- 
III (2006-2008) and SK-IV (2008-2012). The previously re- 
ported Super-K result fl8l . which was consistent with the 
expectation of tau neutrino appearance from oscillation, was 
based on the data from SK-I alone. Since that time the anal- 



ysis has been improved to increase its sensitivity and the data 
set used has been expanded to also include SK-II and SK-III, 
thereby almost doubling its size. As the total data set covers 
the period between 1996 and 2008, it comprises 2806 days of 
live-time. 

In order to predict the rate of both the tau signal and at- 
mospheric background, a full Monte Carlo (MC) simulation 
is used both to predict the neutrino interactions inside the de- 
tector and to model the response of Super-K itself. Three- 
dimensional neutrino fluxes for Vj, and v e produced in atmo- 
spheric showers are taken from the flux calculation of Honda 
et al. Il22l . The fluxes are oscillated with a custom code E3l 
which takes into account all relevant path lengths, energies 
and matter effects using our current knowledge of three-flavor 
neutrino oscillation parameters. The oscillation parameters 
used are EH E3: Amj 2 = 2.1 x 10~ 3 eV 2 , tsm\ x = 7.6 x 
10~ 5 eV 2 , sin 2 2023 = 10, sin 2 20 12 = 0.84, <5 CP = 0. For the 
values of 613 the recent Daya Bay ifTBI and RENO 1 16 1 results 
are combined in a weighted average and we use sin 2 20b = 
0.099. The interactions of the v^, V e , and oscillation-produced 
V T s with the nuclei of water molecules inside of the Super-K 
detector are modeled with the NEUT [2 26 1 neutrino inter- 
action code. Finally, a GEANT3 [27 1 based detector Monte 
Carlo is used to simulate Super-K itself. More detailed de- 
scriptions of this software can be found in |f2)- 

For the purposes of this analysis it is important to under- 
stand some details of the neutrino interaction model. The 
NEUT code models the known neutrino-nucleon interac- 
tions including quasi-elastic scattering, single meson produc- 
tion, coherent pion production, and deep-inelastic scattering 
(DIS). All and v e interactions are simulated. Addition- 
ally, charged-current (CC) V T interactions are simulated and 
added to the sample using weighting based on the oscillation 
probabilities. Neutral current (NC) interactions are assumed 
to be unaffected by oscillations. Tau leptons are decayed us- 
ing TAUOLA (Version 2.6) |28|- Since the distribution of de- 
cay particles depends on the polarization of the tau lepton, 
a polarization model [29 1 is incorporated into NEUT. At the 
relevant neutrino energies selected by this analysis, the CC in- 
teractions contain a high percentage of DIS (46%) with the v T 
signal interactions alone being 56% DIS . In the calculation 
of the cross sections of DIS, the GRV94 1 30 1 parton distribu- 
tion functions are used, and additional corrections to make the 
DIS cross-sections match smoothly with the resonance region 
as developed by Bodek and Yang are also applied [31 1. More 
details of the DIS implementation can be found in [2|. 

The signature of oscillation-induced tau neutrinos in the at- 
mospheric flux is the detection of the decay of tau leptons in 
the Super-K detector. As the leptonic decays of the tau look 
on the whole very similar to normal CC DIS interactions from 
Vjj and v e , an analysis is developed which attempts to select 
the hadronic decays. 

In order to select the tau events, we first identify high en- 
ergy events contained in the inner detector by requiring that 
there is no appreciable activity in the outer detector, the in- 
teraction is in the fiducial volume (the distance to the nearest 
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wall > 200 cm), and the event has more than 1.3 GeV of vis- 
ible energy. The selection efficiencies for this set of cuts are 
81% for the v T CC signal and 23% for the background events 
respectively. The presence of the extra pions from hadronic 
tau decay which come from a heavy object can statistically 
separate the signal from the normal and v e CC and NC 
background. 

In order to further separate the signal from the background, 
a set of variables which show differentiation between the two 
samples are used as the inputs to a feed-forward neural net- 
work (NN). The NN is configured using the TMVA pack- 
age |[32l with seven inputs nodes, one hidden layer with 10 
nodes and one output node. Exclusive training and testing 
samples are selected from the Monte Carlo sets to avoid bias 
and test for overtraining. 

The variables used are (a) the log base ten of the total visible 
energy of the event, (b) the particle ID of the maximum energy 
ring in the event, (c) the number of decay electron candidates 
in the event, (d) the maximum distance between the primary 
interaction and any decay-electron found from a pion or muon 
decay, (e) the clustered sphericity of the event in the center of 
mass system, (f) the number of possible Cherenkov ring frag- 
ments, and finally (g) the fraction of total number of photo- 
electrons in the events carried by the first ring. The agreement 
between downward going data and Monte Carlo (where no tau 
signal is expected) along with the signal expectation is shown 
in Fig. [1] 

All of the oscillation-induced tau neutrinos will come from 
below due to differing path lengths in the earth. In order avoid 
encoding such up-down biases into the expectations, and to 
select events based solely on their topology, the training is per- 
formed by weighting the oscillation probabilities of all events 
based on their energies alone, not their direction. In this way, 
all oscillation probabilities are correct on average, but upward- 
going and downward-going events are treated the same in the 
training process. This technique has the added benefit of not 
setting the weights of the down-going signal events to zero, 
thus preserving Monte Carlo statistics. The training is per- 
formed such that an NN output which is near 1 .0 signifies that 
the event is tau-like, while events near 0.0 are non-tau like. 
After training, the NN is found to efficiently separate the tau 
appearance signal from the background of other atmospheric 
neutrino interactions. 

When acting on the events passing the pre-selection cuts, 
75% of the signal events (60% of the total) and only 26% of 
the background events (6% of the total) remain when events 
with a NN output of greater than 0.5 are considered. In this 
"tau-like" sample, NC background makes up 26% of the sam- 
ple and is an important remaining background. Table [[[further 
displays the fractional breakdown of the interaction modes in 
the sample. In order to extract maximum information from the 
event samples, instead of cutting on the NN output, the output 
of the NN is combined with the zenith direction of the event 
into a probability distribution functions (PDF) and is used to 
jointly fit the tau and background components. 

An example of the two-dimensional distributions of the NN 
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FIG. 1. The distributions of the seven input variables to the NN. The 
downward going data (points) are overlaid with the downward going 
atmospheric MC (solid line). Also shown is the tau signal Monte 
Carlo (shaded). The content of each plot is explained in the text. The 
tau normalization is arbitrary. 



Interaction Mode 

CC v e 

CC Vfi 

CCv t 

NC 



NN < 0.5 
781.4 (0.40) 
1070.2 (0.55) 
12.4 (0.01) 
95.2 (0.05) 



NN > 0.5 All 

381.3 (0.46) 1162.7 (0.42) 

200.2 (0.24) 1270.4 (0.46) 
37.2 (0.04) 49.7 (0.02) 

209.3 (0.25) 304.4(0.11) 



TABLE I. The fractional breakdown of interaction modes of both 
the expected signal (CC v T ) and background for the SK-I period. 
For fitting purposes the entire sample is used in the analysis, but the 
NN enhanced (NN >0.5) and depleted (NN<0.5) signal selections are 
shown here to demonstrate the effect of signal and background sep- 
aration. For each sample, the number of selected SK-I MC events 
is shown scaled by the 1489 days of SK-I live-time. The fractional 
breakdown by interaction mode of each sample is shown in paren- 
theses. 



output versus the direction of the detected events used to dis- 
criminate signal from background is shown in Fig. [2] Distribu- 
tions for both oscillation-generated taus on the left, and other 
atmospheric background on the right are shown. The vertical 
axes of these two-dimensional distributions contain the output 
of the NN and reflects how tau-like the event is (NN output 
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FIG. 2. Histograms of the PDFs of both tau signal (left) and atmospheric background (right). The vertical axis is the output of the NN, the 
horizontal axis the cosine of the event zenith direction. Upward going events are to the left, downward going to the right. The tau signal 
appears in the upward-going tau-like region. 



near 1.0 tau-like, 0.0 non tau-like). The horizontal-axis is the 
cosine of the reconstructed zenith angle of the event which is 
determined by an energy-weighted sum of the ring directions 
in the event. The tau events (left panel) are indicated as tau- 
like by the NN and come from below (cos(0) near -1.0) as ex- 
pected. Atmospheric neutrinos (right panel) on the other hand 
are primarily non-tau-like and come from both above and be- 
low. In fact, it can be seen that these events are depleted in the 
upward-going direction due to their oscillation into (mostly 
non-interacting) tau-neutrinos. By varying the relative nor- 
malization of the two distributions both the amount of tau ap- 
pearance and the overall background level can be adjusted. 

PDFs for each run period for both signal and background 
are built out of two-dimensional histograms prepared from the 
Monte Carlo, with the probability density following the nor- 
malized bin contents. Then, an un-binned likelihood fit of the 
data is done to the sum of the signal and background PDFs 
varying the normalization between them. It is necessary to 
perform an un-binned fit as statistics of bins in the full two di- 
mensional space would be quite low. The result of the fit is a 
normalization factor on the signal and the background which 
tell us how many tau interaction are needed to be consistent 
with our data set. Separate PDFs are produced for SK-I, II and 
III, and each data set is fit to its appropriate MC set. The data 
sets are fit both individually for each run period, and jointly 
together. 

Although the technique employed here is more sophisti- 
cated than that of ifTSl . it is also more sensitive to some sys- 
tematic errors since large numbers of background events re- 
main in the non-signal regions of the fitting space which were 
previously removed by cuts. By training the NN to recog- 
nize tau interactions, the NN also learns to effectively separate 
quasi-elastic (QE) from multi-pi and DIS interactions. This is 
because the DIS events tend to have many extra pions in them, 
and thus look more like the tau signal. The DIS portion of the 



interactions thus form a large part of the background in the 
signal region and we therefore explicitly take into account un- 
certainties in the DIS normalization in the fit. 

The average neutrino energy in the DIS interactions in our 
sample is 14 GeV and the cross section is not known to bet- 
ter than the 10% level at that energy. We also know that the 
application of the Bodek-Yang corrections BP tends to sup- 
press our DIS interactions at higher energies by about 5%. For 
this reason, the DIS error is introduced into the fit as a 10% 
gaussian error constraint. After the fit is completed it is found 
that the amount of DIS is increased from its nominal value 
by 10% at the best fit point. If the fit is performed with no 
constraint on the DIS fraction at all, then the DIS fraction fits 
14% higher than the nominal value. 

The fit is performed on each data period separately, and is 
also performed jointly with all data periods being fit at the 
same time. In the case of finding the exact expectation as pre- 
dicted by the Monte Carlo these factors would be 1.0. When 
the data periods are fit together, the tau normalization is found 
to be 1.42 ±0.35 (stat) with the background normalization 
0.94 ±0.02 (stat). When fit separately, the tau normaliza- 
tions are found to be 1 .27 ± 0.49 (stat), 1 .47 ± 0.62 (stat ) and 
2. 16 ± 0.78 (stat) for SK-I, SK-II and SK-III respectively. 

It is also instructive to examine the results of the combined 
fit graphically. Binned projections of the fitted results can il- 
lustrate the quality and features of the fit. Figure [3] shows the 
projections in zenith for both tau-like (NN output > 0.5) and 
non-tau-like (NN output < 0.5) events, along with the pro- 
jections in NN output for both upward going (cos(0) > 0.1) 
and downward-going (cos(0) < 0.1) events. In these plots the 
PDFs have been rescaled to the fitted normalization values. 
The fitted tau signal is shown in grey. Good agreement is seen 
in all distributions. As expected, tau events are observed as 
an excess of tau-like events in the upward-going direction. In 
these plots the PDFs and data sets from all of the run periods 
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have been combined together. 
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FIG. 3. Fit results showing projections in the NN output and zenith 
angle distribution for tau-like (NN>0.5), upward going (cos(6) > 
0.1), non-tau-like (NN<0.5), and downward going (cos(0) < 0.1) 
events for both the two-dimensional PDFs and data. The PDFs and 
data-sets have been combined from SK-I through SK-III in this fig- 
ure. The fitted tau signal is shown in grey. 

There are 28 uncertainties which are a subset of those used 
in the Super-K three-flavor atmospheric neutrino analysis. A 
more detailed description of them can be found in ll33ll . The 
systematic errors for the analysis are divided into two sets. 
The first set, which describes errors on the tau expectation it- 
self, plays no role in comparing the observed number of events 
with the no-tau appearance hypothesis and does not affect the 
significance of this quoted result. However, this set is used 
to quote an error on the expected number of events and in- 
cludes uncertainties in v T cross section, and any uncertainty 
that would increase both the signal and the background in a 
way that does not change the significance of the reported re- 
sult. Detector biases on selection and fitting are included in 
these uncertainties but are quite small compared to the tau 
cross-section error, the largest being a 5% error on the detector 
energy scale. The error on the tau cross-section was made by a 
comparison of NEUT 11261 with several other models, looking 
in particular at the differences between NEUT and the cross- 
section model by Hagiwara et al ll29l . Another comparison 
between cross-section models was recently completed by the 
authors of [20 1 and gave similar results. 

The second class of errors includes those that would af- 
fect the observed signal but not the background or otherwise 
would cause the significance of the measured normalization 
to change when doing the fit. There are five such errors: the 
up/down flux ratio, the horizontal/vertical flux ratio, the K/n 
flux ratio, the NC/CC cross-section ratio, and the up/down en- 



ergy scale difference in the detector. All of these errors can 
in principle cause a measured change in the significance of 
the signal normalization. In the current analysis the dominant 
error on the signal was the NC/CC ratio changing the fitted 
number of events of about ±7% due to the relatively large 
percentage of NC background in the signal region. 

Also included in the errors which can change the mea- 
sured results and significance are those due to variations in 
the known oscillation parameters. For this study they are var- 
ied within the la limits of the combined SK-I+SK-II+SK-III 
atmospheric oscillation analysis results assuming the normal 
hierarchy 11241 . The Am^ 2 is varied between 1.92 x 10~ 3 and 
2.22 x l(T 3 eV 2 , sin 2 20 2 3 is varied between 0.93 and 1.0. The 
013 values are varied within our combined Daya Bay ITT31 and 
RENO [16| results of sin 2 20 13 = 0.099 ± .014. The use of 
non-zero 0n results in a 13% reduction of the fitted normal- 
ization as three flavor oscillations produce high energy up- 
ward going electron neutrinos which add to the upward-going 
background, thus decreasing the needed number of tau neutri- 
nos to explain the signal region. However, the variation in 613 
around this central value results in less than a 1 % change in 
the fit result. For this analysis, we set the value of 8cp to zero. 
Varying the value of 8cp results in, at most, a 1.3% difference 
in the number of fitted taus, and we neglect this uncertainty. 
The systematic errors are summarized in Table [ill 



Systematics Uncertainties for Vi 


normalization 


+ % 


- % 


Super-K atmospheric v oscillation errors 






28 error terms 


(expected events) 


13.4 


14.7 


5 error terms 


(observed events) 


7.9 


8.5 


Tau neutrino cross section 


(expected events) 


25.0 


25.0 


Oscillation parameters 


(observed events) 


5.4 


1.3 



TABLE II. Summary of systematic uncertainties for both the ex- 
pected and observed number of v T events. The errors for each cate- 
gory including that of the oscillation parameters have been added in 
quadrature. 



Including and combining the observed (+9.6 -8.6% ) 
and expected (+28.4 -30.0%) systematic uncertainties sep- 
arately, the fitted value of the tau normalization is 1 .42 ± 
0.35 (stat) +q|j (sys). After rescaling the Monte Carlo by all 
fitting factors and correcting for efficiency, the observed num- 
ber of fitted events over the entire running period is calculated 
to be 180.1 ±44.3 (stat) +^-| (sys) events. This is to be com- 
pared to an expectation of 120.2 (sys) interactions in the 
fiducial volume if no fitting factors are applied. 

The observed number of events is converted to the signifi- 
cance level at which we can reject the no-tau-appearance hy- 
pothesis. The measured signal normalization (1.42 and its as- 
sociated statistical and systematic errors) is compared with the 
case of no v T appearance, which would have a normalization 
of zero. An asymmetric gaussian centered at 1 .42 is prepared 
and the integral of the PDF below zero is calculated. The p- 
value is 6.2 x 10~ 5 which corresponds to a significance level 
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of 3.8<7. A significance of 2.7(7 is expected for the nominal 
expected signal. The larger measured significance is a con- 
sequence of the fact that more signal was measured than ex- 
pected. The DIS fraction is fit with a 10% increase over its 
nominal value, and is correlated with the tau normalization. 
Because of this, not only is the fitted tau normalization lower 
than it would be without this error, but the error on the tau 
normalization is larger than it would otherwise be due to the 
presence of the correlated DIS error, thus slightly reducing the 
measured significance. It should be noted that if the inverted 
hierarchy is chosen instead of the normal one when calculat- 
ing the oscillation probabilities, the expected number of On- 
induced upward going electrons is reduced, approximately in 
half, resulting in a somewhat higher fitted value (1.56) and a 
correspondingly higher significance. 

In summary, we find that the Super-Kami okande atmo- 
spheric neutrino data is best described by neutrino oscilla- 
tions that includes tau neutrino appearance in addition to the 
overwhelming signature of muon neutrino disappearance. By 
a neural network analysis on the zenith angle distribution of 
multi-GeV contained events, we have demonstrated this at a 
significance of 3.8(7. 
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